/*
 *---------------------------------------------------
 * Random distributions from Numerical recipes in C
 * 1989
 *
 * Typed in by Marcelo Walter
 *---------------------------------------------------
*/

/****************************************************

	INCLUDES

****************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/****************************************************

	DEFINES

****************************************************/
#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 134456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349
#define MINE_PI 3.141592654

/****************************************************

	PROTOTYPES	

****************************************************/
float gammln();
float gasdev();

/**************************************************

        FUNCTIONS

****************************************************/
/*
 * Returns a uniform deviate between 0.0 and 1.0
 * Set idum to any negative value to initialize
 * or reinitialize the sequence
 * pg. 210 of the book.
*/
float ran1(int *idum)
{
  static long ix1,ix2,ix3;
  static float r[98];
  float temp;
  static int iff=0;
  int j;

  if(*idum < 0 || iff == 0){
    iff=1;
    ix1=(IC1-(*idum))%M1;
    ix1=(IA1*ix1+IC1)%M1;
    ix2=ix1%M2;
    ix1=(IA1*ix1+IC1)%M1;
    ix3=ix1%M3;
    for(j=1;j<=97;j++){
      ix1=(IA1*ix1+IC1)%M1;
      ix2=(IA2*ix2+IC2)%M2;
      r[j]=(ix1+ix2*RM2)*RM1;
    }
    *idum=1;
  }
  ix1=(IA1*ix1+IC1)%M1;
  ix2=(IA2*ix2+IC2)%M2;
  ix3=(IA3*ix3+IC3)%M3;
  j=1 + ((97*ix3)/M3);
  if(j>97||j<1){
    printf("Error in ran1.This should not happen!\n");
    exit(0);
  }
  temp=r[j];
  r[j]=(ix1+ix2*RM2)*RM1;
  return (temp);
}


/*
 *---------------------------------------------------------------
 * Returns as a floating point number an integer value that is
 * a random deviate drawn from a Poisson distribution of mean
 * "xm", using ran1(idum)
 * as a source of uniform random deviates.
 * Pg. 294 from the book
 *---------------------------------------------------------------
*/
float poidev(float xm,int *idum)
{
  /* oldm is a flag for whether "xm" has changed since last call */
  static float sq, alxm, g, oldm=(-1);
  float em,t,y;
  float gammln();
  
  if(xm < 12.0){
    if(xm != oldm){
      oldm = xm;
      g = exp(-xm);
    }
    em = -1;
    t=1.0;
    do{
      ++em;
      t *= ran1(idum);
    } while (t>g);
  } else{
    if (xm != oldm){
      oldm = xm;
      sq = sqrt(2.0*xm);
      alxm = log(xm);
      g = xm*alxm - gammln(xm+1.0);
    }
    do{
      do{
	y = tan(MINE_PI*ran1(idum));
	em = sq * y + xm;
      } while( em < 0.0);
      em = floor(em);
      t = 0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
    } while( ran1(idum) > t);
  }
  return em;
 }


/*
 *------------------------------------------------------------
 * Returns a normally distributed deviate with zero mean and
 * unit variance, using ran1(idum) as the source of uniform
 * deviates
 * pg. 217 from the book
 *------------------------------------------------------------
*/
float gasdev(idum)
int *idum;
{
  static int iset=0;
  static float gset;
  float fac,r,v1,v2;
  float ran1();
  
  if(iset==0){
    do{
      v1=2.0*ran1(idum)-1.0;
      v2=2.0*ran1(idum)-1.0;
      r=v1*v1+v2*v2;
    }while (r>=1.0);
    fac=sqrt(-2.0*log(r)/r);
    
    gset=v1*fac;
    iset=1;
    return v2*fac;
  }
  else{
    iset=0;
    return gset;
  }
}
	

/*
 * Natural log of the gamma function
 * Pg 168 from the book
 *
 */
float gammln(xx)
float xx;
{

  double x, tmp, ser;
  static double cof[6]={76.18009173, -86.50532033,
			    24.01409822, -1.231739516,
			    0.120858003e-2, -0.536382e-5};
  int j;

  x = xx-1.0;
  tmp = x + 5.5;
  tmp -= (x+0.5)*log(tmp);
  ser = 1.0;
  for(j=0; j<=5; j++){
    x += 1.0;
    ser += cof[j]/x;
  }	
  return -tmp + log(2.50662827465*ser);
}













